Novel mitochondrial gene rearrangements pattern in the millipede Polydesmus sp. GZCS‐2019 and phylogenetic analysis of the Myriapoda

Abstract The subphylum Myriapoda included four extant classes (Chilopoda, Symphyla, Diplopoda, and Pauropoda). Due to the limitation of taxon sampling, the phylogenetic relationships within Myriapoda remained contentious, especially for Diplopoda. Herein, we determined the complete mitochondrial genome of Polydesmus sp. GZCS‐2019 (Myriapoda: Polydesmida) and the mitochondrial genomes are circular molecules of 15,036 bp, with all genes encoded on + strand. The A+T content is 66.1%, making the chain asymmetric, and exhibits negative AT‐skew (−0.236). Several genes rearrangements were detected and we propose a new rearrangement model: “TD (N\R) L + C” based on the genome‐scale duplication + (non‐random/random) loss + recombination. Phylogenetic analyses demonstrated that Chilopoda and Symphyla both were monophyletic group, whereas Pauropoda was embedded in Diplopoda to form the Dignatha. Divergence time showed the first split of Myriapoda occurred between the Chilopoda and other classes (Wenlock period of Silurian). We combine phylogenetic analysis, divergence time, and gene arrangement to yield valuable insights into the evolutionary history and classification relationship of Myriapoda and these results support a monophyletic Progoneata and the relationship (Chilopoda + (Symphyla + (Diplopoda + Pauropoda))) within myriapod. Our results help to better explain the gene rearrangement events of the invertebrate mitogenome and lay the foundation for further phylogenetic study of Myriapoda.

In recent decades, with the development of molecular biology, a relatively new field of molecular analysis based on mitochondrial and transcriptome datais flourishing. In contrast to this traditional morphology view, several molecular studies contradicted the Dignatha clade and supported Symphyla + Pauropoda group formed Edafopoda (Andreas et al., 2012;Dong et al., 2012;Fernandez et al., 2018;Gai et al., 2006;Rehm et al., 2014). Therefore, the relationship among the four classes of myriapod is still controversial and the major source of conflict is between molecular and morphological phylogeny.
The millipedes (Diplopoda) is an important component of the modern terrestrial ecosystem due to its important role in the decomposition of organic matter. Hitherto, the Diplopoda contained more than 18,000 species worldwide, which are distributed in most parts of China (Golovatch & Liu, 2020;Jiang & Chen, 2018).
Although Diplopoda is the third most diverse class of Myriapoda, there is no widely accepted consensus about the classification and phylogenetic relationship. With the development of molecular biology technology, a new era of phylogenetic analysis of phylogeny has been opened in the early 1990s, and a large number of analyses of millipedes have been published (Brewer et al., 2013;Dong et al., 2016;Lavrov et al., 2002;Liu et al., 2017;Means et al., 2021;Qu et al., 2020;Wesener et al., 2010;Zhao et al., 2020). However, due to the limitation of taxon sampling and the lack of mitochondrial genome data, the previous phylogenetic studies failed to solve the relationship of millipedes.
The mitochondrial genomes of metazoans exhibit variation in many characteristics, such as length, tRNA secondary structure, gene rearrangement, and structure of control regions (Boore, 1999;Mukundan et al., 2020;Shen et al., 2017Shen et al., , 2020. Studying the variation in these characteristics can discover the evolutionary relationship between taxa with a high and/or low classification level. Among them, gene arrangements are relatively complex and diverse, which can become a source of information for system evolution analysis. Furthermore, it also affects the process of mRNA transcription, substitution, and processing. In recent years, the mitochondrial genome rearrangement has been widely studied focusing on phylogenetic relationship and rearrangement mechanism (Feng et al., 2021;Gong et al., 2020;Li et al., 2019Li et al., , 2020Powell et al., 2020;Tyagi et al., 2020;Wang et al., 2020;Zhang et al., 2020). In addition, the high rearrangement rate makes the Myriapoda an ideal group to study the interaction between gene rearrangement and phylogenetic relationship. For example, the studies discussed the gene arrangement of Myriapoda on phylogenetic inference but did not prove the universality of this mechanism in the same order species (Gai et al., 2008;Lavrov et al., 2002). Some studies found that the gene arrangement pattern was a sound molecular evidence supporting the Helminthomorpha clade (Brewer et al., 2013;Dong et al., 2012), but they did not elaborate the evolutionary implications of gene arrangements in the Myriapoda. Several common models have been used to explain the gene rearrangement events in the current animal mtDNA, for example: recombination models involved in DNA strand breaks and recombination (Lunt & Hyman, 1997); the Tandem duplication-random loss (TDRL) model is commonly used to support gene tandem replication and random loss (Moritz & Brown, 1987); and the TDNL model supports gene tandem replication and non-random loss (Lavrov et al., 2002). However, the gene rearrangement phenomenon may not be explained by one of the above-mentioned mechanisms alone for some species. Therefore, it is necessary to conduct comparative evolutionary studies on mitogenome rearrangements to accurately identify the mechanisms leading to the rearrangements.
In the present study, we sequenced the complete mitochondrial genome of a millipede, Polydesmus sp. GZCS-2019 (P. GZCS-2019) (Diplopoda: Polydesmidae), and described the genome-scale gene rearrangement events of the mitogenome, providing independent molecular evidence to explore the phylogenetic relationship of Myriapoda. To build a better phylogenetic relationship and understand the evolutionary significance of gene arrangement in Myriapoda, the other 27 complete mitochondrial genomes of the Myriapoda (8 from Chilopoda, 13 from Diplopoda, 2 from Symphyla, and 1 from Pauropoda) and 3 outgroup species were used in this study. Meanwhile, we combine phylogenetic analysis, divergence time, and gene arrangement to yield valuable insights into the evolutionary history and classification relationship of Myriapoda.

| Specimen collection and DNA extraction
Two specimens of P. GZCS-2019 were collected from Chishui of Guizhou Province in China (28°24′25″N, 105°57′17″E) in August 2019. Morphological identification of specimens was mainly referred to as "PICTORIAL KEYS TO SOIL ANIMALS OF China" (Yin, 1998) and all specimens were stored in anhydrous ethanol in the Chongqing Key Laboratory of Animal Biology, Chongqing Normal University. Total DNA was extracted from the dehydrated muscle tissues using the TaKaRa MiniBEST Universal Genomic DNA Extraction Kit Ver.5.0 (TaKaRa Biotech).

| Mitochondrial genome sequencing and assembly
The entire mitogenome of P. GZCS-2019 was sequenced on the Illumina HiSeq TM platform with paired ends of 300-500 bp. The raw paired reads were quality trimmed using FastQC v0.11.4 (www. bioin forma tics.babra ham.ac.uk/proje cts/fastq c/) with default parameters. Finally, yielded 10G raw reads (coverage 3-5×) and clean sequence reads were assembled in the NOVOPlasty (https://github. com/ndier ckx/NOVOP lasty) (Nicolas et al., 2016) using sequences from each of the 23 mitochondrial genes of the closest relative available from NCBI as mapping reference, with the default parameter.

| Sequence analysis and gene annotation
The online tool MITOS (http://mitos2.bioinf.uni-leipz ig.de/index.py) was used to perform gene annotation, and the annotation results were verified by the BLAST program from the NCBI website (https:// blast.ncbi.nlm.nih.gov/Blast.cgi) (Donath et al., 2013). Then, the abnormal start codon and stop codon were determined based on the comparison with other millipedes. The relative synonymous codon usage (RSCU) was obtained using MEGA 7.0 (Kumar et al., 2015), which was calculated using PCG with incomplete codons removed.
The ribosomal RNA genes were determined according to the location of adjacent tRNA genes and comparison with other Myriapoda mitogenomes from NCBI. The strand asymmetry was calculated using the following formula: AT-skew = (A − T)/(a + T); GC-skew = (G − C)/ (G + C) (Perna & Kocher, 1995). The online mitochondrial visualization tool Organellar Genome DRAW (Marc et al., 2013) was used to draw a graphical map of the mitochondrial genome. The secondary cloverleaf structure and the locations of tRNAs were examined with tRNAscan-SE 1.21 (Lowe & Chan, 2016).

| Phylogenetic reconstruction
The mitochondrial genomes used for phylogenetic analysis in this study were all from GenBank, including 24 species of Myriapoda and 3 species of outgroup (1 Decapoda species and 2 Hexapoda species).
The species information is shown in Table 1. This phylogenetic analysis is based on 37 genes, including 13 protein coding genes (PCG), 2 ribosomal RNA genes (rRNAs) and 22 transfer RNA genes (tRNAs).
The sequences above were aligned by ClustalW method in MEGA 7 (Kumar et al., 2015), with the default parameters. The Gblocks version 0.91b (Castresana, 2000) with the default parameters setting was used for filtering of poorly aligned regions. The aligned sequences of each gene were concatenated using Sequence Matrix v1.7 (Castresana, 2000).
Phylogenetic trees were constructed using the following three datasets: (1) 13 PCGs matrices consisting of 9140 nt; (2) 13 PCG and 2 rRNA matrices composed of 10,884 nt; (3) 13 PCGs, 2 rRNAs, and 10 tRNA matrices composed of 11,459 nt. For these three datasets, the best fitting model GTR + I + G was selected by jModelTest 2 (Darriba et al., 2012) for maximum likelihood (ML) and Bayesian inference (BI) analysis. The ML analysis was assembled in PhyML 3.0 (Stéphane & Olivier, 2003) with fast likelihood-based method and performed 1000 repetitions. Bayesian analyses were carried out using MrBayes 3.1.2 (Ronquist & Huelsenbeck, 2003) under the best-fit models with 10,000,000 generations in two runs of eight chains each and each one was sampled every 200 generations with a burn-in of 25%. Trees inferred prior to stationarity were discarded as burn-in, and those remaining were used to construct a 50% majority rule consensus tree. All phylogenetic trees were viewed and edited using Figtree v1.3.1 (http://tree.bio.ed.ac.uk/softw are/figtree).

| Divergence time estimation
Beast v1.8.4 (Drummond et al., 2012) was used to estimate the divergence time, using the Bayesian analysis method. At the same time, Beauti v1.8.3 was used to generate the beast XML file using uncorrelated lognormal distribution relaxed clock model and the Yule speciation was used to prior process the tree. Two fossil constraints were used in this study: the oldest uncontested terrestrial animal Pneumodesmus newmani (421-426 Mya) (Shear et al., 1998) and the oldest terrestrial myriapod body fossil Rhyniella praecursor (407-411 Mya) (Wilson & Anderson, 2004). The GTR + I + G model was used to estimate time, and after a burn-in of the initial 50% cycles, divergence times were sampled once every 1000 generations from 100 million Markov Chain Monte Carlo (MCMC) iterations.
The treeAnnotator v1.6.1 (BEAST software) was used to annotate the sampled trees, and the Figtree v1.3.1 was used to conduct the visualization. The ESSs were used to determining the Bayesian statistical significance of each parameter in TRACER v1.5 (ESS > 200) (Rambaut & Drummond, 2003).

| Genome structure, organization, and composition
The complete mitogenome sequence of P. GZCS-2019 is a closed circular molecule with a size of 15,036 bp ( Figure 1 and Table 1). In addition, the gene content also conforms to the typical characteristics of other Diplopoda species, including 13 PCGs (cox1-3, nad1-6, nad4L, cob, atp6, and atp8), 2 rRNA genes (rrnS and rrnL), 22 tRNA genes, and a control region, and all genes are encoded on the heavy (+) chain ( Figure 1 and Table 2). Moreover, the mitogenome contains 351 bp intergenic spacer sequences, distributed in 19 regions, ranging in size from 1 to 174 bp (Table 2), and there is a 27 bp overlap between genes in five locations, showing five pairs of overlapping genes: atp8/atp6, rrnS/trnV, trnP/nad4L, nad4L/nad4, and trnH/ nad5, of which the longest 7 bp overlap is located between trnL1 and rrnL, nad4L and nad4. Furthermore, the whole mitogenome of TA B L E 1 Summary of mitogenomic sequence information used in the present study

| PCGS and codon usage
In the mitogenome of P. GZCS-2019, the size of the PCGs region is 9882 bp, and all the PCGs genes are encoded on the + strand (  Figure S4).

| Skewness, transfer RNAs, and ribosomal RNAs
The nucleotide composition of the mitogenome of P. GZCS-2019 is as follow: A (25.3%), T (40.9%), G (24.2%), and C (9.7%) ( Table 3).  except trnS1 and trnS2 without dihydrouridine (DHU) arm ( Figure   S2). In general, such deletion of DHU arm in the secondary structure of trnS1 and trnS2 was considered a common condition in the Diplopoda mitogenome (Brewer et al., 2013;Dong et al., 2016  The rrnS (1054 bp) gene located between Control region and trnV, and the rrnL (808 bp) gene located between trnV and trnL1 are encoded on + strand (Table 1 and Figure 1). The A+T content (70.2%) of the rRNA genes is higher than the whole genome (66.1%) (Table 3), with a negative AT-skew −0.056 (Table 2) and whose structural diagram is shown in Figure S3.

| Control regions
The largest non-coding region of the mitochondrial genome is usually presumed as the control region and is heavily biased toward A+T which were identified as initiation sites for replication and transcription (Boore, 1999;Cameron et al., 2007;Jeffrey, 1999;Shadel & Clayton, 1997;Taanman, 1999;Wei et al., 2013). However, the poly A-stretches at the 5′ and 3′ end of the ancestor L. polyphemus are not observed in the other three Polydesmida species ( Figure S6).
Therefore, we speculate that this event is responsible for the reverse of strand transcription direction observed in this Polydesmida order and further experiments are needed to clarify this speculation.
However, further experiments are needed to clarify this speculation.
Tandem duplication-random loss mechanism was widely used to explain the translocation of mitochondrial genes, the trnT translocation phenomenon in this study can be explained by this theory, that occurring in the region between nad4L and trnP, followed by deletions of redundant genes resulting in trnT-nad4-nad4L (Figure 2b2).
In contrast, the inversion of trnC-trnQ referred to the transversion of trnI-trnQ, which was more in line with the recombination model Furthermore, we also found the inversion of the entire side of a genome trnP,trnQ,trnC,and trnY) and the translocation of trnT and the inversion of trnC-trnQ could be proposed as common events about gene order in Polydesmida lineage ( Figure S5). The duplication-nonrandom loss was also detected in the Symphyla species (Gai et al., 2006(Gai et al., , 2008, which reinforce the sister relationship with Diplopoda. These results of the regular gene arrangement in Myriapoda provide useful information for the phylogenetic inference of advanced groups.

| Divergence time
Understanding the origin and evolutionary history of myriapods is crucial for interpreting the colonization and evolution of arthropods on land. Hitherto, Myriapoda is inferred to have colonized land in the Early Cambrian, substantially predating body or trace fossil evidence (Giribet & Edgecombe, 2019;Lozano-Fernandez et al., 2016;Wilson & Anderson, 2004). In this study, the Bayesian divergence times showed that the splitting of the ancestral lineages of the Progoneata and Chilopoda from a common ancestor occurred during the Wenlock period of Silurian, slightly earlier than the oldest millipedes and centipedes fossil records in the Silurian, which was similar to many previous studies (Rosa et al., 2016;Wilson & Anderson, 2004), suggesting that both had experienced the same period of relative stasis period as the plants prior to this (Giribet & Edgecombe, 2019;Lozano-Fernandez et al., 2016;Minter et al., 2017). Then, the Progoneata clade split into Dignatha (Diplopoda + Pauropoda) and Symphyla during the early Silurian to lower Devonian. The results showed that the first split of myriapod occurred between the subclass Chilopoda and other subclasses, rather than between Symphyla and other subclasses, which align with the morphology classification results that supports the monophyly of Progoneata ( Figure 5).
In Chilopoda, the split time of Pleurostigmophora and Notostigmophora was from the early Silurian to middle Triassic, slightly earlier than the oldest fossil chilopods in the Late Silurian (Wilson & Anderson, 2004), which was consistent with some previous studies and these were representatives of the Chilopoda (Bonato et al., 2015;Chipman et al., 2014;Giribet & Ed Gecombe, 2013). Moreover, this study also concluded that the split time of Scolopendromorpha and (Lithobiomorpha + Geophilomorpha) was from middle Devonian to early Jurassic and the divergence time of Lithobiomorpha and Geophilomorpha was from early Carboniferous to early Cretaceous ( Figure 5).
In Diplopoda, the divergence time of the two infra-  Figure 5).

| CON CLUS ION
In this paper, we report the complete mitogenome of P. GZCS-2019 (Diplopoda: Polydesmidae) with a novel genome-scale rearrangement phenomenon. We deduce the genome-scale F I G U R E 5 The divergence time estimation of the major myriapod lineages using the Bayesian relaxed molecular clock method in BEAST from two fossil constraint ages based on the best scoring maximum-likelihood tree. Node bars indicate 95% CIs of the divergence time estimate "duplication + (non-random/random) loss + recombination (TD (N\R) L + RC)" model resulted in a novel mechanism of gene rearrangement for the published Polydesmida mitogenome. The deletion of the DHU arm of trnS1 and trnS2 was considered a common condition in the Polydesmida mitogenome. The phylogenetic analysis supported the monophyletic of Diplopoda, providing evidence for the higher-level relationships within it. Meanwhile, we combine phylogenetic analysis and divergence time to yield valuable insights into the evolutionary history and classification relationship of Myriapoda and these results support a monophyletic Progoneata and the relationship (Chilopoda + (Symphyla + (Diplopoda + Pauropoda))).
Since the mitochondrial gene rearrangement events in Myriapoda contain genetic information related to the phylogenetic evolution of species, it is necessary to conduct in-depth research and use the genetic information revealed by gene rearrangement to better solve these controversial phylogenetic problems. However, due to the lack of taxon samples, there are still many limitations in this study.
Therefore, it is necessary to collect more species of Myriapoda for more in-depth and systematic research.

ACK N OWLED G M ENTS
The authors sincerely thank all the crews for their help with the manuscript writing and data analysis. This work was supported by the

Science and Technology Research Program of Chongqing Municipal
Education Commission (Grant No. KJQN201900502).

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
DNA sequences: GenBank (Accession Numbers MZ677220).